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Analysis  of  intensity  structure  of  the  ISI 
model  in  the  FAST2D  hydrocode 

Andrew  J.  Schmitt 

Laser  Plasma  Branch,  Plasma  Physics  Division 
Naval  Research  Lab  Wash.  DC  20375 

Consider  a  lens  illuminated  by  an  ISI  beam 
and  focused  onto  a  target.  Fourier  optics  tells  us 
that  the  spatial  distribution  of  the  electric  field  at 
the  lens  is  a  direct  analog  of  the  wavenumber  spec¬ 
trum  of  the  electric  field  at  the  target.  For  the  mo¬ 
ment,  assume  that  the  beam  is  normal  to  the  target. 
The  extent  of  the  field  in  wavenumber  space  at  the 
target  is  due  entirely  to  the  geometry  of  the  focus¬ 
ing  optics:  the  diameter  of  the  lens  D,  the  focal 
length  of  the  lens  /,  or  the  F/#  of  the  lens  F=flD.  As¬ 
sume  further  that  the  intensity  profile  of  the  beam  at 
the  lens  is  such  that  it  uniformly  covers  the  lens  sur¬ 
face,  is  uniform  in  a  statistical-time-average  sense, 
and  is  characterized  by  the  number  of  spatial  co¬ 
herence  zones,  Nz,  that  cover  the  lens  surface.  The 
resulting  envelope  of  the  electric  field  or  intensity 
at  the  target  will  have  a  width  Rspot  =  NzFAq,  or  Nz 
times  the  diffraction  limited  spot  size  of  the  entire 
lens,  F\ o-  The  spatial  width  of  the  coherence  zone 
at  the  lens  surface  is  D/Nz.  This  electric  field  struc¬ 
ture  at  the  lens  results  in  a  spectral  width  of  the 
electric  field  coherence  zones  at  the  target  surface 
of  Akfzln/FAo  as  well  as  the  total  spectrum  width 
there,  kM=NzAkj  =2%NzIFXq.  Thus,  the  spectrum  va¬ 
ries  from  -kM/2<k<+kM/2  where  k  is  the  wave- 
number  of  the  field  at  the  target. 

(The  preceding  limits  assume  that  the  angle 
0  subtended  by  the  lens  is  small  enough  that 
sin0~0.  More  generally,  the  spatial  coherence  zone 
at  position  -D/2<x<D/2  and  width  d  at  the  lens  will 
result  in  an  electric  field  at  the  target  at  wave- 
number  k=kosin(Bx)  (where  0x=atan(x#)  and  kQ  is 
the  laser  wavenumber  2k/'Kq)  and  width 
Aki=ko(sin[atm((x-d/2)/f)]-sin[atan((x+d/2)/f)]).) 

Multiple  Beams: 

We  can  consider  multiple  beams  in  the  limit 
that  they  are  immediately  adjacent  to  one  another; 
this  is  equivalent  to  shortening  the  F/#  of  the  incid¬ 
ent  lens  optic  to  Fmui^flNbD=FINb  where  F  is  the 


LENS  PLANE  ELECTRIC  FIELD  (SPACE) 


TARGET  PLANE  ELECTRIC  FIELD  (WAVENUMBER) 


Fig.  1.  Focusing  optics  produce  the  Fourier  trans¬ 
form  of  the  incident  beam  on  the  target. 

F /#  of  the  single  lens.  This  results  in  the  width  of 
the  ^-spectrum  of  the  electric  field  at  the  target  in¬ 
creasing  by  the  factor  Nb>  while  the  incoherent  zone 
width  Akj  remains  the  same  (the  total  number  of  co¬ 
herence  zones  is  multiplied  by  Nb).  This  leads  to 
higher  frequency  structure  in  both  the  electric  field 
and  the  intensity.  When  this  shorter  wavelength 
structure  is  averaged  over  a  region  larger  than  the 
smallest  wavelengths,  smoothing  of  the  resulting 
beams  occurs.  We  will  consider  this  effect  in  more 
detail  later. 
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1.1-1  D:  Intensity  structure: 

Given  the  ^-spectrum  of  the  electric  field  at 
the  target,  the  intensity  can  be  found  as  its  convolu¬ 
tion: 

I(x)  =  E\x)E{x)  =  J  dke~ik  xE*(k  )J  dk  eikxE(k ) 

=  \dk\  dk  eikxE*  (k  )E(k  +  k  ) 

(1) 

In  the  simulation  code,  the  integral  is  effectively  re¬ 
placed  by  a  summation  (the  Discrete  Fourier  Trans¬ 
form)  as  in: 


=  O) 

j 

where  £  =  k?JAk{  varies  from  -Nz+1  to  Nz-1,  and  the 
limits  of  the  sum  over  j  are  from  j=\-Nz/2- 
MIN(f,0)  to  j=Nz/ 2-MAX(^,0). 

Now  if  Ej  are  complex  random  values  r\  + 
z>2  where  r\  &  r 2  are  instances  of  a  Gaussian  ran¬ 
dom  distribution  of  width  unity,  then  Ej  is  a  vector 
in  the  complex  plane  with  amplitude  a  and  a  uni¬ 
formly  distributed  random  phase  between  0  and  2%. 
For  any  £^0,  E*Ej+i  will  also  have  a  random  phase 
and  amplitude  <&> 2.  Summing  n  of  these  numbers 
will  (on  average)  give  a  vector  of  amplitude  nl/2<a>2 
with  a  random  phase.  For  each  £,  the  number  of 
sums  given  by  Eqn.  (2)  is  simply  Nz- \£\.  Thus  the 


wavenumber  index  (k/Ak) 

Fig.  2:  The  solid  dark  line  shows  the  time  averaged 
simulation  results  after  1000  coherence  times,  and 
the  broken  red  line  shows  the  analytic  predictions, 
for  Nz=75,  Iavg=l  (eqn.  4).  The  wavenumber  index 
corresponds  to  l  in  the  text. 


(average)  amplitude  of  the  intensity  from  (2)  can  be 
written: 


(mi) = llE*Ej+\ = <«>2 


(3) 


\J  / 

where  <>  denotes  an  average  element  from  the  sta¬ 
tistical  ensemble.  A  special  case  exists  for  £-0  since 
the  products  E*Ej  are  all  real  and  positive;  the  sum 
there  is  Nz<&2>.  Finally,  since  the  £=0  term  also 
gives  the  average  value  of  the  intensity  Iavg,  we 
have:  Iavg  =  \dx<l(x)>/\dx  =  <I(k=0)>=Nz<a2>;  thus 
<a2>=Iavg/Nz.  However,  note  the  subtle  distinction 
between  <a>2  and  <&2>  —  they  are  not  equal  as  the 
operations  of  ensemble-averaging  and  raising  to  an 
power  do  not  commute.  For  the  probability  distri¬ 
bution  function  we  consider  here,  p(I)  =  exp(- 
///avg),  which  gives  </>=/avg  and  thus  <a2>=Iavg/Nz  as 
noted  above;  but  as  p(\E\)  =  2(|£'|/4vg)exp(-|£,|2//avg), 
<|£|>  =  <d>  =  V(7t<a2>)/2  =  (l/2)(7t/avg/Az)i/2  (see  Ap¬ 
pendix  A  and  fig.  2) 


(4) 


where  1 <=k/kM.  As  we  can  see  from  fig.  2,  this  for¬ 
mula  matches  the  simulation  model  very  well,  when 
averaged  over  a  suitable  time,  and  when  the  simula¬ 
tion  is  set  up  to  resolve  all  the  modes  in  the  intensi¬ 
ty  spectrum. 

Before  leaving  this  section,  I  want  to  point 
out  two  things.  First,  note  that  the  effect  of  adding 
multiple  beams  is  to  simply  to  increase  the  number 
of  incoherence  zones  Nz  by  the  factor  Nf,  (number 
of  beams),  thus  reducing  each  Fourier  component 
of  the  intensity  by  the  square  root  of  the  number  of 
beams,  as  expected.  However,  multiple  beams  also 
multiply  the  extent  of  the  intensity  spectrum  by  Nj,. 

Secondly,  note  that  the  actual  (time  average) 
amplitude  of  an  intensity  perturbation  at  a  given 
wavelength  \-=2n/\kf\  is  twice  the  quantity  given  in 
eqn.  4,  because  one  must  add  both  the  positive  and 
negative  wavenumber  contributions. 

I.2A-1D  The  problem  of  incomplete  resolu¬ 
tion:  Inadequate  Fourier  resolution 

It  is  a  common  problem  in  the  simulation 
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that  we  cannot  resolve  all  the  modes  generated  by 
the  laser  beam.  For  instance,  if  a  simulation  has  a 
transverse  size  (X[en)  smaller  than  the  spot  size  of 
the  beam,  the  modes  with  wavelength  X[en<X<Rspot 
(and  others  not  fitting  the  periodicity  of  the  mesh) 
cannot  be  resolved.  In  general,  in  a  simulation  of 
size  X ien  with  Nx  points  along  x  (grid  size 
A x=Xien/Ny)  and  periodic  boundary  conditions 
(needed  by  ISI  and/or  Fourier  transforms,  and  in 
any  case  more  general  than  predefined  amplitude  or 
gradient  boundaries),  the  Fourier-space  grid  will 
have  extent  -kx/2<k<kx/2  and  resolution  Akx  where 
Akx  =  2K/Xien  and  kx=2%/Ax=NxAkx.  This  resolution 
may  not  match  the  resolution  needed  for  the  field, 
e.g.  in  general,  Ak^Akj.  In  particular  we  consider 
here  the  case  of  inadequate  resolution  such  that 
Akx>Akj  (i.e.,  Xien<Rspot )•  Then  each  Fourier  mode 
in  the  code  must  represent  Nr=AkxJAkf=Rspot/Xien 
modes  (assume  for  the  moment  Nr  is  integer;  the  re¬ 
sult  will  generalize  to  the  non-integer  case).  A  typi¬ 
cal  value  of  Nr  in  the  code  would  be  7.5 
(#*/>0f=75Op.m  [Nike],  X/e„=100pm). 

There  are  many  ways  to  model  this.  One  is 
to  simply  (linearly)  interpolate  the  real  spectrum  at 
the  places  where  we  do  resolve  the  field  (i.e., 
kja=iAkx,  i=-Nx,...,Nx).  Another  way  is  to  use  the 
field  averaged  over  the  cell  A£xin  Fourier  space. 

Consider  the  first  case,  that  of  simple  inter¬ 
polation  (i.e.,  the  amplitude  of  the  electric  field  is 
“sampled”  at  the  point  in  fc-space  at  which  it  is 
centered).  Denoting  E  as  the  field  resolvable  in  the 
simulation,  we  have:  t{kxi)=E{ki=k^).  The  change 
in  the  number  of  modes,  however,  implies  that  the 
average  intensity  would  change,  unless  we  renor¬ 
malize  the  field  amplitudes,  because  we  no  longer 
have  Nz  modes  but  instead  Nz/Nr  modes.  In  this  case 
<E>i=<a’>2=(nIavglX)N/Nz  and  eqn  (3)  becomes: 

('(*»•» = -jNjNr-mf  - <5> 

Consider  the  second  way,  and  average  the 
field  amplitude  over  the  Fourier  grid  cell  A  kx: 
E(kxi)=\&kx  E(kx  )/ Akx=(Akj/ Ak^)XE(kxi+nAki), 
where  the  summation  is  over  the  Nr  nearest  modes. 
Assuming  that  the  sum  is  again  over  statistically 
identical  but  mutually  incoherent  fields, 
<£(kxi)>~Nrm<E(kx)>=Nrv2d.  &  is  still  determined  by 
the  average  intensity  and  the  (total)  number  of 


modes:  d=(nIavg/4Nz)U2.  However,  we  have  also  re¬ 
duced  the  number  of  terms  in  the  sum  of  eqn  (2)  to 
Nz'  =  Nz/Nr.  Thus  Eqn.  (3)  becomes: 


=  l>K(<S2) 


l-k 


Nz/Nr 


(6) 


...which  is  identical  to  eqn.  (5).  So  in  each  case  the 
necessity  of  normalizing  the  electric  field  amplitude 
to  produce  the  correct  average  intensity,  in  tandem 
with  our  insistence  that  each  electric  field  compon¬ 
ent  is  statistically  identical,  leads  us  to  the  same  so¬ 
lution.  Compared  to  the  perfectly  resolved  case 
(eqn  (4)),  the  individual  intensity  fluctuations  are 
larger  by  a  factor  of  -Jkf  .  These  formulas  (5)  or 
(6)  represent  what  is  in  the  FAST2D  code  now. 

This  increase  in  mean  intensity  fluctuation 
level  is  the  result  of  preserving  the  total  rms  fluc¬ 
tuation  and  average  intensity  while  keeping  the  an- 
satz  that  each  electric  field  mode  amplitude  is  iden¬ 
tical.  One  way  around  this  increase  in  the  am¬ 
plitude  would  be  to  simply  interpolate  the  mode 
amplitudes,  but  drop  the  constraint  of  producing  the 
correct  average  intensity.  We  would  find  the  intens¬ 
ity  spectrum: 


(wcyt))=pjNr-mi 


l-kl  ^ 


avg 


'NzNr 


(6a) 


The  correct  average  intensity  can  still  be  recovered 
by  simply  adding  the  appropriate  (constant)  amount 
to  the  DC  mode  (k=0):  instead  of  </(&= 0)>= 
Nz<a2>/Nn  add  the  amount  /avg(l-l/Ar).  In  this  meth¬ 
od,  the  individual  intensity  fluctuations  are  smaller 
by  a  factor  of  -Jkf  compared  to  the  resolved  case 
(eqn  (4)),  and  a  factor  of  Nr  smaller  than  the  cur¬ 
rent  methods  (5)  and  (6). 

Finally,  the  alteration  in  individual  mode 
amplitude  produced  by  the  methods  (5 -6a)  can  be 
avoided  by  directly  using  the  intensity  spectrum 
given  by  eqn.  (4)  no  matter  what  resolution  we 
might  have.  Then  (if  needed)  add  the  appropriate 
DC  contribution  to  preserve  the  average  intensity. 
This  method  preserves  the  average  intensity  and  the 
individual  intensity  fluctuation  level,  at  the  expense 
of  reducing  the  total  rms  intensity  fluctuation  level 
(as  there  are  fewer  modes  represented). 

(In  the  last  two  methods,  there  should  prob¬ 
ably  be  a  random,  statistical  nature  to  this  add-on 
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2sin(£A/2)/M  as  a  function  of  kA,  and  the  dotted 
line  is  its  envelope,  2/kA. 

constant  —  since  this  DC  term  is  actually  being 
made  up  of  all  from  those  Nr-l  modes  that  we  failed 
to  resolve  during  the  interpolation.  Also,  it  may  be 
appropriate  to  reduce  the  mean  variance  of  the  indi¬ 
vidual  mode  amplitude  (by  ~Nr-V2)  to  reflect  the  ef¬ 
fect  of  averaging  the  adjacent  unresolved  modes.) 

I.2B-1  D:The  problem  of  incomplete  resolu¬ 
tion:  Inadequate  spatial  resolution 

Another  common  problem,  particularly 
when  using  multiple  beams,  is  not  having  enough 
spatial  resolution  to  resolve  the  high  frequency  spa¬ 
tial  structure  (remember,  the  smallest  intensity 
wavelength  is  expected  to  be  of  order  FXq,  where  F 
in  this  case  is  represents  the  total  beam  cluster).  Let 
us  denote  A  as  our  smallest  resolvable  spatial  extent 
(n.b.:  A  may  be  the  same  size  as  —or  larger  than— 
the  mesh  spacing  Ax).  Then  the  resolved  intensity  at 
any  scale  due  to  this  minimum  resolution  is: 

-  *+A/2  jc+A/2 

4W  =  -  jdx'IE(x)  =  ^Ik  jdx'elkx' 

x-A/2  k  x-A/2  gy 

_  ^^gifcc|2sinM/2| 

where  the  subscript  “R”  refers  to  the  resolved  in¬ 
tensity  and  the  subscript  “E”  refers  to  the  exact  in¬ 
tensity,  and  4  are  the  Fourier  coefficients  of  the 
“exact”  intensity.  Thus,  averaging  the  exact  intensi¬ 
ty  over  the  distance  A  is  equivalent  to  multiplying 


the  Fourier  coefficients  of  the  intensity  by  the  factor 
2sin(M/2)/M.  (see  fig.  3).  Note  that  at  kA=n%  for 
integer  n  greater  than  zero,  the  multiplier  is  zero. 

What  about  the  high-frequency  end  of  the  k 
spectrum  that  can’t  be  resolved  in  the  simulation, 
e.g.  \k\>nJX[en?  These  values  of  k  are  aliased  to  val¬ 
ues  of  k  that  do  appear  in  the  simulation.  The  rule 
for  such  aliasing  is  that  the  actual  k  is  mapped  to  a 
resolvable  k  given  by: 

**es  =  i*m-|(*MOD*m)-|*m|  (8) 
where  MOD  is  the  modulo  function.  Thus  the  in¬ 
tensity  at  each  resolvable  wavenumber,  k,  is  made 
up  of  the  weighted  sum  of  the  intensity  at  wave- 
numbers  k,  ku-k,  kM+k,  2kM-k, ... 

II.A-2D:  The  actual  (two  dimensional) 
illumination  spectrum 

Now  we  concern  ourselves  with  the  actual 
illumination  spectrum  with  two  transverse  dimen¬ 
sions  (i.e.,  the  “real  world”  case).  Again,  we  assume 
that  the  k-spectrum  of  the  electric  field  at  the  target 
is  known  (as  it  is  simply  equivalent  to  the  spatial- 
distribution  of  the  electric  field  at  focusing  lens). 
And  again,  as  in  the  development  of  eqn.  (1),  the 
intensity  can  be  found  as  its  convolution: 

I(x,y)  =  E*(x,y)E(x,y) 

\ldk,dk>ei-x"'l’yBkx.ky)  (9) 

E*  (kx ,  ky  )E(kx  +kx,ky+ky) 
and  now  the  discrete  summation  is  also  in  two 
transverse  directions: 

m,kn)  =  ^  ZElmEj+e,m+n  (10) 

m  j 

This  summation  is  determined  by  the  number  of 
modes  in  the  sum,  since  each  summand  is  statisti¬ 
cally  identical.  The  number  of  sums  is  given 
(geometrically)  by  the  area  of  overlap  of  two  circles 
of  width  knf/2,  separated  by  a  distance 
I k|=V(k;c2+ky2).  This  area  of  overlap  is  the  number 
of  sums  in  eqn.  (10)  (in  units  of  A k2  (we  assume 
Ak=Akx=Aky)),  and  is  given  by: 
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(11) 


where  mk/kM  varies  from  0  to  1.  So  if  each  indi¬ 
vidual  electric  field  amplitude  is  a,  the  amplitude  of 
the  intensity  fluctuation  at  a  given  K  - 
flkx2+ky2)/kMis: 

{Kk'.kj)  = 

(13) 

Again,  the  special  case  exists  at  K=0,  where  the  sum 
is  coherent  and  the  mean  intensity  is  the  average  in¬ 
tensity:  Ns(K^0)<a2>=(k^f/Ak)2K<a2>/4=InVg.  This  de¬ 
termines  <a2>:  <d2>=4IavsjnNz2,  (Nz=kM/Ak)  and  by 
extension  (see  discussion  previous  to  eqn.  (4)  and 
Appendix  A)  <a>2=Iavg/Nz2.  So  the  average  intensity 
fluctuation  at  wavelength  K  is  given  by: 

(/(«:))  =  2sin_1  [k] - 2rcJ] L  —  K*2  (14) 

where  Nz=kM/Ak  is  the  number  of  coherence  zones 
measured  across  a  ID  slice  of  the  lens  (the  total 
number  of  coherence  zones  is  ~NZ2).  Note  that  for 
small  (but  finite)  K,  </(  k=$Q)>=>  -Jtc/{2Nz  ) .  This 


fig.  4:  The  2D  spectrum  of  intensity  fluctuations  as 
a  function  of  wavenumber  K-sqrt( k^ + ky^)/kM  (in 
red).  The  ID  spectrum  (eqn.  4)  is  shown  for  com¬ 
parison;  both  use  Nz=l  for  comparison  purposes. 


spectrum  as  a  function  of  K  is  very  similar  to  the  ID 
case,  except  that  it’s  smaller  (per  mode)  by  a  factor 
~Nz1/2 

II.B-2D:  Properly  Accounting  for  2D  effects 
with  a  1 D  treatment 

The  major  difference  between  the  2D  and  the 
ID  treatment  of  ISI  is  that  the  number  of  modes 
with  large  \k\  is  larger  in  2D:  the  number  of  modes 
with  a  given  wavenumber  |fc|  increases  linearly  as 
~\k\.  One  (arbitrary)  way  to  take  this  into  account  in 
ID  is  to  preserve  the  relative  <7nns  (compared  to  the 
full  2D  distribution)  at  a  given  wavenumber.  The 
standard  deviation  of  the  2D  distribution  is  calculat¬ 
ed  straightforwardly  (using  eqn.  (14),  and  assuming 
a  round  focusing  lens): 


The  integrand  is  thus  simply  proportional  to  the 
arms  for  each  mode  in  the  2D  distribution.  For  com¬ 
parison,  the  ID  distribution  gives  (using  eqn.  4): 

/  £t\2  X  '  /  Tf1.\\2 


w =  £<w 


jdk(I(k))2+ jdk(I(k))2  jtk 


=  2  $dk(I(k))2  /&k 


l^(i-^)  =  ^42vg(i-i^+^- 

V*. 


7?  J2 

16  1avg 
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(/Of.))- 


(18) 


fig.  5:  The  function  converting  ID  imprint  to  2D,  as 
a  function  of  wavenumber  K(c.f.  eqn.  17). 


The  total  Gnns  in  2D  is  the  same  as  in  ID,  but  since 
there  are  many  more  modes  in  2D,  the  typical  Grms 
in  any  given  individual  mode  is  much  lower. 

The  ID  distribution  can  be  “renormalized”  to  2D 
by  requiring  that  (La)  the  individual  ID  mean-in¬ 
tensity  fluctuation  in  each  mode  is  equal  to  the  2D 
mean-intensity-fluctuation;  and  ( i.b )  the  ID  total 
orms  is  the  same  as  the  2D  Grms.  Alternatively,  we 
can  (ii)  create  a  new  ID  distribution  by  averaging 
the  2D  distribution  in  k-space  to  obtain  a 
“representative”  distribution. 

Method  (i)  would  be  used  to  scale  a  ID 
(transverse)  simulation  to  a  2D  (transverse)  simula¬ 
tion;  in  this  case  the  imprint  at  each  wavenumber  of 
a  ID  (transverse)  simulation  would  be  scaled  by  the 
multiplier  composed  of  eqn(14)/eqn(4): 


WO>2d  =  r 


2 


/r - 2 sin  I[(c]-2yVl-y2 


ID 


(17) 

(see  figure  5);  this  multiplier  is  of  order  unity  ex¬ 
cept  for  the  Nzl/2  factor.  Using  this  expression  to 
renormalize  the  mode  amplitudes  and  summing 
them  over  the  resulting  2D  distribution  also  yields 
the  correct  total  Grms. 

Method  (ii),  creating  a  new  ID  distribution  by 
averaging  the  2D  distribution  in  |k|-space,  could  be 
used  if  we  wished  to  generate  a  2-D  like  modal 
Grms  in  a  ID  (transverse)  simulation  directly.  Re¬ 
quiring  that  the  intensity  fluctuation  amplitude  at  k* 
reflect  the  total  2D  fluctuation  level  at  k,  the  ID 
“equivalent”  of  the  2D  field  would  be: 


/r-2sin 


-l 


kx-2kx 


This  has  been  normalized  by  requiring  that  the  total 
Grms  be  the  same  as  in  2D.  It’s  not  clear  how  to  get 
this  spectral  shape  simply  by  generating  the  elec¬ 
tric  fields  as  we  did  before  —  the  electric  field  com¬ 
ponents  are  no  longer  independent  of  kx.  Probably 
the  easiest  way  is  to  assemble  the  field  with  the  typ¬ 
ical  ID  intensity  distribution  (sec.  I),  then  apply  a 
lD-to-2D  transform  in  Fourier  space  by  using  the 
ratio  of  eqn  (18)  to  eqn  (4)  (see  also  fig.  6): 


{^Kx^2D-mock  _  1 

[4**1 

/r-2sin  kx-2kx-s) 

M) 

k{\  -  k ) 

(19) 

(which  is  just  ^tcNzkx  times  eqn.  (17)).  As  before, 
missing  modes  (due  to  incomplete  resolution  in 
Fourier  space)  can  be  accounted  for  by  normalizing 
the  resulting  field  to  ensure  that  the  total  G^  is 


fig.  6:  The  normalization  function  (eqn  (19))  ap¬ 
plied  to  a  ID  distribution  spectrum,  in  order  to  get 
a  2D  intensity  distribution  spectrum. 


preserved. 
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Appendix  A 

The  distribution  function  of  the  electric  field  magni¬ 
tude  can  be  derived  from  the  distribution  function 
of  the  intensity  using  the  cumulative  distribution 
function  (e.g.,  B.W.  Lindgren  Statistical  Theory 
(2nd  Ed.)  p.  69  (1968),  or  any  basic  text  on  statis¬ 
tics).  In  general,  we  can  equate  the  two  cumulative 
distribution  functions  (e.g.,  F(X)),  which  are  de¬ 
fined  from  the  (density)  distribution  functions  (e.g., 
p(x))  as: 

x 

F(X)=  J dxp(x )  (A.l) 

— oo 

In  the  case  under  consideration  here,  we  know  the 
intensity  distribution  function  p(I),  and  we  wish  to 
find  the  electric  field  distribution  function  Jp(|£'|). 
Since  p{I)  =  exp[-///av5]//avg,  eqn.  (A.l)  gives 
F](J)=  1-exp  [-I/Iavg]/Iavg>  which  is  equal  to  the  cu¬ 
mulative  electric-field  distribution:  F|£|(|£|)=l- 


exp[-\E\VIaVg]/Iavg.  (Iavg=<F.)  The  probability  densi¬ 
ty  function  p(|£|)  is  found  from  the  cumulative  dis¬ 
tribution  through  differentiation  as: 

=Sk=2_ (A.2) 

^1^1  -yfiavg 
(This  has  been  normalized  using  the  requirement 
that  the  total  probability  is  unity).  Furthermore,  the 
average  value  of  |2J|,  <|£)>  is  found  as  the  first  mo¬ 
ment  of  the  probability  density  function: 

oo  / 

(A-3) 

o 

We  note  that,  although  the  average  value  of  the  in¬ 
tensity  distribution  is  (by  definition)  Iavg,  the  aver¬ 
age  value  of  the  electric  field  distribution  is  not 
simply  /avg1/2,  but  is  instead  only  about  89%  of  that 
value. 


Attachment  A 

The  IDL  function  implementing  the  ISI  model,  similar  to  that  used  in  the  FAST2D  code,  and  used  for  the 
comparisons  above  (fig.  2). 


function  ISIS2,  Nbeams,  Nmodes,  NMX,  seed,  AVINT=AVINT ,  debug=debug,  $ 
ntl=ntl,  verbose= verbose,  eef=eef 


;  NAME: 

;  isis2 

y  CALLING  SEQUENCE: 

;  intensity  =  isis2 (Nbeams ,  Nmodes,  Nmax,  SEED) 

;  PURPOSE: 

;  generate  an  ISI  profile  of  length  NMX,  return  in  array  IFlux 
\  FUNCTION  RETURNS: 

;  IFLUX:  array  of  size  NMX  that  is  output  as  isi  intensity 


?  PARAMETERS: 
;  Nbeams : 

7  NMODES : 


NMX: 

seed: 


Number  of  beams  that  we'll  be  adding  (?) 

Integer  number  of  ISI  modes  (in  k  space)  of  EACH  beam 

curiously,  they  will  be  arranged  from  k=[0  :  NMODES-1] , 

e.g.  biased  to  k.ge.O  in  Fourier  space.  Not  that  this  should  matter. 

Integer  size  of  output  Iflux  vector 

Random  number  seed,  (should  be  provided  so  that  it  isn’t  reset  to  same 
value  on  every  call) 


7  OPTIONAL  PARAMETERS 

7  Avint:  Average  intensity  (Def  =  1.0) 

;  ntl :  Number  of  total  modes  per  resolved  mode  (def  =  1.0) 

;  (In  practice,  this  will  be  Y__focal /Y_length) 

7  verbose:  Verbose  flag  (def=0) ?  if  nonzero,  debugging  notices  get  printed 

7  debug:  debug  flag  (def  =  0) 

7  eef:  The  electric  field 

HISTORY: 

;  programmed  circa  04  Nov  1999  :  in  preparation  for  APS/DPP  99  confernce 
7  modified  from  isis.pro  on  04  Jan  2000  :  testing  out  new  ideas  of 
7  multiple  beams  &  averaging. 

if  n_elements (avint)  eg  0  then  avint  =  1.0 

if  n_elements (debug)  eq  0  then  debug  =  0 

if  n_elements (ntl)  eq  0  then  ntl  =  1.0 

if  n_elements (verbose)  eq  0  then  verbose  =  0 
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NmodeTot  =  Nmodes  *  Nbeams 

;  find  Np  =  2**MM  such  that  Np  > ( 2 *NmodeTot )  .and.  Np  >=  Nmx 
?  then  Nspace  =  Np/Nmx 

MM  =  FIX (  ALOG (NmodeTot) /ALOG (2. )  +  2  )  >  FIX(  ALOG (Nmx) /ALOG (2 .) ) 
if  (NP  It  Nmx)  then  NP  =  NP*2 

NSPACE  =  NP/NMX  ;  I  think  this  is  guarranteed  to  be  integer 
if  verbose  ne  0  then  begin 

print, 'In  isis2 :  using  vector  of  length  ’,np, 1  to  generate  length  1  ,nmx 
print, '  —  nspace  =  ', nspace, '  averaging  points' 
endif 


construct  a  random  amplitude  distribution  in  Fourier  space 
for  wavenumbers  corresponding  to  indices  NP/4  to  NP/4+NmodeTot 


eef  =  complexarr (np) 

eef(O)  =  CoMPLeX( randomn (seed,  Nmodetot) ,  randomn ( seed,  Nmodetot) ) 

;  add  in  effects  of  other  modes  not  resolved  here 
if  (ntl  gt  1.0)  then  begin 
for  ijk  =  2,  ntl  do  begin 

eef(0)  =  eef(0:*)  +  CoMPLeX (randomn (seed,  Nmodetot),  randomn(seed,  Nmodetot)) 
endfor 

if  ntl  mod  1.0  gt  l.E-4  then  begin 
eef(0)  =  eef(0:*)  +  $ 

(ntl- fix ( ntl )) *CoMPLeX{ randomn (seed,  Nmodetot),  randomn (seed,  Nmodetot)) 

endif 

endif 

eef  =  eef /sqrt (2 . *Nmodetot*ntl) 


;  NOTE:  Normalization  of  l./NP  is  applied  in  directional 
;  (from  IDL  online  help: 

7  A  normalization  factor  of  1/N,  where  N  is  the 

;  number  of  points,  is  applied  during  the  forward 

;  transform. ) 

sav  =  abs (  NP*fft(eef,  -1)  ) rt2 


;  average  every  Nspace  points,  put  it  into  IFlux 
indx  as  indgen(Nmx)  *nspace 
If lux  =  sav (indx)  ;  >0. 

if  Nspace  gt  1  then  begin 

for  ijk  =  1,  Nspace-1  do  begin 

I flux  =  If lux  +  sav (indx  +  ijk) 
endfor 

Iflux  =  Iflux  /  float (Nspace) 
endif 

/  — 

Iflux  as  iflux  *  avint 

if  debug  ne  0  then  stop 

return,  Iflux 

END 
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